function results = fading_mimo_capacity(varargin)
% 衰落MIMO信道容量分析
% 输入参数:
%   'SNR_dB' - SNR范围 (dB), 默认: -5:5:30
%   'antenna_configs' - 天线配置数组, 默认: [1, 2, 4, 8, 16]
%   'num_realizations' - 蒙特卡洛仿真次数, 默认: 500
% 输出:
%   results - 包含容量结果、分布、中断容量等信息的结构体

% 解析输入参数
p = inputParser;
addParameter(p, 'SNR_dB', -5:5:30);
addParameter(p, 'antenna_configs', [1, 2, 4, 8, 16]);
addParameter(p, 'num_realizations', 500);
parse(p, varargin{:});

SNR_dB = p.Results.SNR_dB;
antenna_configs = p.Results.antenna_configs;
num_realizations = p.Results.num_realizations;
SNR_linear = 10.^(SNR_dB/10);

% 添加路径
addpath('../Common');

% 获取颜色定义
colors = color_definitions();

fprintf('=== 衰落MIMO信道容量分析 ===\n');

% 瑞利衰落信道容量
fading_capacity_csi = zeros(length(SNR_dB), length(antenna_configs));
fading_capacity_no_csi = zeros(length(SNR_dB), length(antenna_configs));

for ant_idx = 1:length(antenna_configs)
    Nt = antenna_configs(ant_idx);
    Nr = antenna_configs(ant_idx);
    
    fprintf('仿真 %dx%d 瑞利衰落MIMO...\n', Nt, Nr);
    
    for snr_idx = 1:length(SNR_dB)
        snr = SNR_linear(snr_idx);
        
        % 蒙特卡洛仿真
        capacities_csi = zeros(num_realizations, 1);
        capacities_no_csi = zeros(num_realizations, 1);
        
        for real = 1:num_realizations
            % 生成瑞利衰落信道
            H_fading = sqrt(0.5) * (randn(Nr, Nt) + 1i * randn(Nr, Nt));
            
            % 已知CSI容量
            [U, S, V] = svd(H_fading);
            singular_values = diag(S);
            capacities_csi(real) = water_filling_capacity(singular_values, snr);
            
            % 未知CSI容量
            capacities_no_csi(real) = log2(det(eye(Nr) + (snr/Nt) * H_fading * H_fading'));
        end
        
        % 遍历容量 (平均值)
        fading_capacity_csi(snr_idx, ant_idx) = mean(capacities_csi);
        fading_capacity_no_csi(snr_idx, ant_idx) = mean(capacities_no_csi);
    end
end

% 绘制衰落信道容量
figure('Name', '衰落MIMO信道容量', 'Position', [150, 150, 1200, 800]);

subplot(2,2,1);
for ant_idx = 1:length(antenna_configs)
    plot(SNR_dB, fading_capacity_csi(:, ant_idx), ...
         ['-', colors(ant_idx)], 'LineWidth', 2);
    hold on;
end
grid on;
xlabel('SNR (dB)');
ylabel('容量 (bps/Hz)');
title('衰落MIMO容量 (已知CSI)');
legend(arrayfun(@(x) sprintf('%dx%d', x, x), antenna_configs, 'UniformOutput', false));

% CSI知识的影响
subplot(2,2,2);
ant_idx = 3; % 4x4 MIMO
plot(SNR_dB, fading_capacity_csi(:, ant_idx), 'b-', 'LineWidth', 2);
hold on;
plot(SNR_dB, fading_capacity_no_csi(:, ant_idx), 'r--', 'LineWidth', 2);
grid on;
xlabel('SNR (dB)');
ylabel('容量 (bps/Hz)');
title(sprintf('%dx%d MIMO: CSI影响', antenna_configs(ant_idx), antenna_configs(ant_idx)));
legend('已知CSI', '未知CSI');

% 容量分布分析 (特定SNR)
snr_fixed = 15; % dB
snr_idx = find(SNR_dB == snr_fixed, 1);

% 生成更多样本来分析分布
num_samples = 10000;
capacity_samples = zeros(num_samples, 1);
Nt = 4; Nr = 4;

for i = 1:num_samples
    H_sample = sqrt(0.5) * (randn(Nr, Nt) + 1i * randn(Nr, Nt));
    capacity_samples(i) = log2(det(eye(Nr) + (10^(snr_fixed/10)/Nt) * H_sample * H_sample'));
end

subplot(2,2,3);
histogram(capacity_samples, 50, 'Normalization', 'pdf');
grid on;
xlabel('容量 (bps/Hz)');
ylabel('概率密度');
title(sprintf('容量分布 (SNR=%d dB, %dx%d)', snr_fixed, Nt, Nr));

% 中断容量分析
outage_probabilities = [0.01, 0.05, 0.1, 0.2];
outage_capacities = zeros(length(outage_probabilities), 1);

for out_idx = 1:length(outage_probabilities)
    sorted_capacities = sort(capacity_samples);
    outage_idx = round(outage_probabilities(out_idx) * num_samples);
    outage_capacities(out_idx) = sorted_capacities(outage_idx);
end

subplot(2,2,4);
plot(outage_probabilities, outage_capacities, 'o-', 'LineWidth', 2);
grid on;
xlabel('中断概率');
ylabel('中断容量 (bps/Hz)');
title(sprintf('中断容量分析 (%dx%d MIMO)', Nt, Nr));

% 组织结果
results.SNR_dB = SNR_dB;
results.antenna_configs = antenna_configs;
results.fading_capacity_csi = fading_capacity_csi;
results.fading_capacity_no_csi = fading_capacity_no_csi;
results.capacity_samples = capacity_samples;
results.outage_probabilities = outage_probabilities;
results.outage_capacities = outage_capacities;

end

function capacity = water_filling_capacity(singular_values, snr)
    % 注水算法容量计算
    num_modes = length(singular_values);
    total_power = length(singular_values); % 归一化总功率
    noise_power = 1/snr;
    
    % 注水功率分配
    water_level = water_filling_mimo(singular_values.^2, total_power, noise_power);
    
    % 计算容量
    capacity = 0;
    for i = 1:num_modes
        power_alloc = max(0, water_level - noise_power/singular_values(i)^2);
        if power_alloc > 0
            capacity = capacity + log2(1 + snr * power_alloc * singular_values(i)^2);
        end
    end
end

function water_level = water_filling_mimo(eigenvalues, total_power, noise_power)
    % MIMO注水算法
    num_modes = length(eigenvalues);
    
    % 排序特征值
    [sorted_eig, sort_idx] = sort(eigenvalues, 'descend');
    
    % 寻找最优水位线
    water_level = 0;
    
    for k = 1:num_modes
        temp_level = (total_power + noise_power * sum(1./sorted_eig(1:k))) / k;
        temp_powers = max(0, temp_level - noise_power./sorted_eig(1:k));
        
        if all(temp_powers >= 0)
            water_level = temp_level;
        else
            break;
        end
    end
end